Genetic diversity and population dynamic of Ziziphus jujuba var. spinosa (Bunge) Hu ex H. F. Chow in Central China

Abstract Phylogeographic research concerning Central China has been rarely conducted. Population genetic and phylogeography of Ziziphus jujuba var. spinosa (also called sour jujube) were investigated to improve our understanding of plant phylogeographic patterns in Central China. Single‐copy nuclear gene markers and complete chloroplast genome data were applied to 328 individuals collected from 21 natural populations of sour jujube in China. Nucleotide variation of sour jujube was relatively high (π = 0.00720, θ w = 0.00925), which resulted from the mating system and complex population dynamics. Analysis of molecular variation analysis revealed that most of the total variation was attributed to variation within populations, and a high level of genetic differentiation among populations was detected (F st = 0.197). Relatively low long‐distance dispersal capability and vitality of pollen contributed to high genetic differentiation among populations. Differences in the environmental conditions and long distance among populations further restricted gene flow. Structure clustering analysis uncovered intraspecific divergence between central and marginal populations. Migrate analysis found a high level of gene flow between these two intraspecific groups. Bayesian skyline plot detected population expansion of these two intraspecific groups. Network and phylogeny analysis of chloroplast haplotypes also found intraspecific divergence, and the divergence time was estimated to occur at about 55.86 Ma. Haplotype native to the Loess Plateau was more ancient, and multiple glacial refugia of sour jujube were found to locate at the Loess Plateau, areas adjacent to the Qinling Mountains and Tianmu Mountains. Species distribution model analysis found a typical contraction‐expansion model corresponding to the Quaternary climatic oscillations. In the future, the distribution of sour jujube may shift to high‐latitude areas. This study provides new insights for phylogeographic research of temperate plant species distributed in Central China and sets a solid foundation for the application of the scientific management strategy of Z. jujuba var. spinosa.

plast genome data were applied to 328 individuals collected from 21 natural populations of sour jujube in China. Nucleotide variation of sour jujube was relatively high (π = 0.00720, θ w = 0.00925), which resulted from the mating system and complex population dynamics. Analysis of molecular variation analysis revealed that most of the total variation was attributed to variation within populations, and a high level of genetic differentiation among populations was detected (F st = 0.197). Relatively low long-distance dispersal capability and vitality of pollen contributed to high genetic differentiation among populations. Differences in the environmental conditions and long distance among populations further restricted gene flow. Structure clustering analysis uncovered intraspecific divergence between central and marginal populations. Migrate analysis found a high level of gene flow between these two intraspecific groups. Bayesian skyline plot detected population expansion of these two intraspecific groups. Network and phylogeny analysis of chloroplast haplotypes also found intraspecific divergence, and the divergence time was estimated to occur at about 55.86 Ma. Haplotype native to the Loess Plateau was more ancient, and multiple glacial refugia of sour jujube were found to locate at the Loess Plateau, areas adjacent to the Qinling Mountains and Tianmu Mountains. Species distribution model analysis found a typical contraction-expansion model corresponding to the Quaternary climatic oscillations. In the future, the distribution of sour jujube may shift to highlatitude areas. This study provides new insights for phylogeographic research of temperate plant species distributed in Central China and sets a solid foundation for the application of the scientific management strategy of Z. jujuba var. spinosa.

K E Y W O R D S
chloroplast genome, genetic structure, nucleotide variation, population dynamic, single-copy nuclear gene markers, Ziziphus jujuba var. spinosa

| INTRODUC TI ON
The evolutionary history of plant species has emerged as a complex interaction of biogeography, climate change, and human forces (Feng et al., 2018;Paola et al., 2017). Phylogeography seeks to understand the comprehensive evolutionary history and distribution of organisms and has become a focus of evolutionary biology (Abbott & Comes, 2004;Riddle, 2010). Of particular importance to Chinese phylogeographic research, the Qinghai-Tibetan Plateau (QTP; Gao et al., 2019;Khan et al., 2018;Qiu et al., 2011), Southwestern China Hou et al., 2018;Zheng et al., 2017), and Northern China (Bai et al., 2010;Zeng et al., 2016) have attracted more attention resulting from the substantial diversity of plant species, the topographically heterogeneous terrain and the complex climate conditions. These studies shed light on the interactions of geography, ecology, and climate in shaping the distribution and genetic pattern of plant species. However, phylogeographic research concerning Central China (such as the Loess Plateau and adjacent areas) has been rarely conducted . Some endemic species are narrowly and concentrated distributed in Central China, such as Xanthoceras sorbifolia Bunge (Zhu, 2016) and Elaeagnus mollis Diels (Du et al., 2020). Furthermore, the Quaternary glaciation cycles had little effect on this area, glaciers only formed in scattered areas of the Loess Plateau and lasted for a relatively short period (Zheng et al., 1998). Phylogeographic study of E. mollis, an endangered deciduous shrub species restricted distributed in the Loess Plateau, revealed multiple potential glacial refugia and allopatric divergence in Central China (Du et al., 2020).
Ziziphus jujuba var. spinosa, also known as sour jujube, is a deciduous shrub plant species belonging to Ziziphus, Rhamnaceae. Z. jujuba var. spinosa is widely distributed in the temperate region of China, especially in areas from the Loess Plateau to the Taihang Mountains, which corresponds to its climate and habitat optimum Zhao et al., 2021). It harbors significant ecological values and is typically used for soil and water conservation in Northern China as the high tolerance to drought and salt stress . The nutritional value of the fruit and the medicinal importance of the seed of sour jujube has kept its economic vitality in China for more than 2000 years (Qu, 1982). It is concluded that Chinese jujube (Z. jujuba Mill) is originated from sour jujube, and the evolutionary path may involve several different patterns (Liu, 1993;Peng, 1991). Genome-resequencing SNP data have been utilized to investigate the genetic structure of sour jujube with Chinese jujube, and the clear separation of sour jujube and Chinese jujube, and the genome evolution and domestication history of sour jujube were highlighted (Guo et al., 2021;Huang et al., 2016;Shen et al., 2021).
As the primitive ancestor of Z. jujuba, sour jujube is considered as a valuable gene pool for the genetic improvement of Chinese jujube . Accompanied by the growth of the plantation of Chinese jujube and the overharvesting of the seeds, the distribution range of sour jujube is undergoing severe fragmentation and decrease (Gao et al., 2008). Previous research on sour jujube mainly focused on taxonomic classification, nutritional and medicinal ingredients, and responses to various abiotic stress (Kang et al., 2008;Li et al., 2017;Liu, 1993;Yan et al., 2017). The genetic diversity and genetic structure of sour jujube have not been well characterized, which hinders the comprehensive understanding and utilization of this important genetic resource. Zhang et al. (2015) investigated the genetic diversity and population structure of sour jujube using SSRs and found high level of genetic diversity (H E = 0.659 and H S = 0.674) and moderate differentiation among populations (F st = 0.091, R st = .068). However, the phylogeographic pattern and population demography of sour jujube corresponding to various climate conditions were still yet to be clarified. Phylogeography study of sour jujube may add acknowledgement of the species' evolutionary history in Central China and provide some clues for the potential migration between this area and other adjacent areas.
In the present study, single-copy nuclear gene markers were utilized to investigate the genetic diversity and population structure of sour jujube distributed across China. Furthermore, chloroplast genome data was used to detect the phylogeographic history of sour jujube. The population dynamic of sour jujube responding to the Quaternary and future climatic changes was also documented.
The aims of this study were as follows: (i) determine the genetic diversity and population structure of Z. jujuba var. spinosa, (ii) describe the phylogeographic history and the potential location of glacial refugia, and (iii) investigate the population dynamic of sour jujube in China. The present study will improve our knowledge of plant phylogeography in Central China and set a solid foundation for future management and utilization of this important germplasm resource.

| Population sampling and DNA extraction
In this study, 328 individuals were sampled from 21 sour jujube populations distributed in China. The geographic information of the populations sampled in this study was illustrated in Figure 1 and Table S1. We also sampled one Chinese jujube (Z. jujuba "Hupingzao") and one Z. mauritiana individual as outgroup in the following analysis. Individuals sampled in the same population were at least 100 m apart. Total DNA was extracted for all sampled individuals from silica gel-dried leaves using the modified cetyltrimethyl ammonium bromide (CTAB) method.

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology; Botany

| Locus selection, amplification and sequencing
Single-copy nuclear gene markers have been widely used in plant population genetic studies, resulting from their cross-utility, high variable, and easily processed features (Du et al., 2015;Hou et al., 2018;Liu et al., 2016;Wang et al., 2019). In this study, we selected five single-copy nuclear gene markers developed in our previous research for amplification and sequencing of all the samples (Hu, Du, Wang, & Han, 2021). Detailed information about the locus was listed in Table 1.
Polymerase chain reaction (PCR) amplification followed the protocols established in . Direct sequencing was performed on an ABI 3730XL DNA analyzer (Applied Biosystems) using the same primers that were used for amplification.

| Genetic diversity and the neutral test
The assembled contigs of each individual were aligned using Clustal X ( Thompson, 1997) and refined manually in Bioedit (Hall, 1999).
The number of segregating sites (S), nucleotide diversity parameters, F I G U R E 1 Sampled populations of sour jujube in this study
mauritiana were used as an outgroup for the MLHKA test.

| Genetic differentiation
The distribution of nucleotide variation among populations at each locus was assessed using AMOVA (analysis of molecular variation) implemented in Arlequin 3.5 (Excoffier et al., 2005) with significance testing evaluated using 1000 permutations of the original data.
Genetic variation was hierarchically partitioned into among populations within species and within populations. Pairwise Wright's fixation index, F st (Wright & Maxson, 1968), was used to measure population differentiation Furthermore, to assign the sampled individuals to different genetic clusters (K), genetic structure analysis was performed in Structure 2.3.4 (Falush et al., 2003) with these five loci. Ten independent runs for each possible value of K from 1 to 10 were performed with a burn-in of 500,000 following 1000,000 MCMC iterations with an admixture model and correlated allele frequencies. The most likely value of K based on the negative natural log-likelihood of the data (LnP(K)) and ΔK was calculated using Structure Harvester (Campana et al., 2011;Evanno et al., 2005). Distruct 1.1 (Rosenberg, 2004) was used to create and visualize the population bar plots.

| Chloroplast genome analysis
To understand the phylogeography history of sour jujube in China, we sequenced and assembled the complete chloroplast genome of 21 individuals (randomly chose one individual from each sampled population). The chloroplast haplotypes were assigned using Dnasp 5.0, and the median-joining (MJ) haplotype network was generated in Popart 1.7 (Leigh & Bryant, 2015)  Four independent MCMC replications were performed, and each replication was run for 100,000,000 steps and sampled every 1000 steps. The first 25% of the sampled trees were discarded as burn-in.

| Coalescent-based population demographic history analysis
To reconstruct the demographic history of sour jujube and the two intraspecific groups revealed in Structure analysis (see below) over time, the historical demography was inferred from Bayesian skyline plot analyses (BSP) implemented in BEAST 1.10 (Suchard et al., 2018). This coalescence-based approach utilized a standard MCMC sampling procedure to evaluate the posterior probability distribution of effective population size backward until the time to the most recent common ancestor of the sampled sequences was reached. All the five loci were used for this analysis. The best-fit substitution model of the concatenated sequence was determined by jModeltest 2 (Posada, 2008). Independent MCMC analysis was run for 1 × 10 8 steps, sampling every 100 steps and discarding 25% of sampled trees as burn-in. For each group, multiple independent analyses were performed with different random seeds to test for convergence, and results of replicate runs were pooled using LogCombiner1.10.4, and skyline plots were visualized with Tracer 1.7.1 (http://tree.bio.ed.ac.uk/softw are/tracer).

| Gene flow analysis
An MCMC maximum likelihood method was utilized to estimate gene flow between two intraspecific groups using Migrate 4.4.3 (Beerli, 2001).
Two important population genetic parameters, θ (four times effective population size multiplied by mutation rate per site per generation) and M (immigration rate divided by the mutation rate) were calculated with the implemented F st estimations. Five independent runs by 10 short chains of 5000 steps and 3 long chains of 50,000 steps were run. Genealogies with a sampling increment of 100 and 10,000 burn-in were recorded.

| Landscape genetics
The biogeographic boundaries were calculated by Monmonier's maxi-

| Impact of environmental factors on genetic structure (isolation by environment)
In order to evaluate the effect of present environmental conditions on the observed pattern of the genetic structure of sour jujube, the correlation between GenD and environmental distance (EnvD) among these sampled populations was calculated by Mantel test using Vegan package implemented in R. Twenty-three environmental variables of the 21 sampled populations were determined for the current climate layers as used in species distribution modeling analysis (see below). EnvD matrix was computed based on the population score on each bioclimatic variable. cvh.ac.cn/) and distribution information listed in Zhang et al. (2015) and Zhao et al. (2021). A total of 519 distribution data were collected, and R package "spThin" was utilized to remove the samples that clustered within 10 km to reduce the sampling deviation impact (Aiello-Lammens et al., 2015), and delete redundant data. To reduce sampling deviation, only one distribution point was kept for each grid (10 × 10 km), and finally 253 effective distribution records were used in the following analysis (Table S3). In Maxent analysis, 80%

| Species distribution modeling
of the distribution localities were used to train the model and 20% were randomly selected to test the model. 10 replicate runs were performed, and replicated run type was set to cross-validate. The random seed was used to ensure consistency in the statistical output between runs. Model performance was evaluated using the area under the receiver operating characteristic curve (AUC). The value of AUC varies between 0.5 and 1, and a higher value means higher model accuracy, although in practice the maximum possible value of AUC is often <1 (Phillips et al., 2006). The Jackknife procedure was

| Nucleotide variation and neutrality test
We successfully sequenced 5 single-copy nuclear gene markers from all the individuals sampled in this study. As shown in Table 2

| Population structure and genetic differentiation
AMOVA analysis was used to investigate the overall distribution of genetic variation (Table 3). Variation among populations ranged from 9.06% at locus SZX2 to 27.96% at locus SZ11 and was significant  further proved the intraspecific divergence within sour jujube. BSP analysis revealed that the effective population size of sour jujube and the two intraspecific groups underwent substantial expansion over time (Figure 3b-d).

| Chloroplast genome analysis
We successfully sequenced and assembled the complete chloroplast genome of 21 sour jujube individuals randomly selected from each sampled population. The length of the complete chloroplast genome varied between 159,399 and 161,279 bp, and GC content ranged from 36.51% to 37.30%. Detailed information about the sour jujube chloroplast genome was listed in Table 5. A total of 21 chloroplast haplotypes (H1-H21) were designated and the outgroup Z. mauritiana contained a private haplotype (H22). The MJ network of the chloroplast haplotypes showed that the chloroplast haplotype native to population YCRC was ancestral, which was also supported by BEAST analysis (Figure 4). Both network and BEAST analysis revealed two intraspecific lineages within sour jujube, which was con-

| Distribution model
Using 10 environmental variables, the population distribution dynamic of sour jujube in different periods in China was modeled with Maxent ( Figure 5). The result showed that the average value of the AUC of training data in different periods was 0.943 while that of the test data were 0.932, meaning higher accuracy of the distribution model (Table S5). The contribution rate of Bio11 was the highest (35.4%), and soil was the least (0.1%) to the model in the Jackknife procedure (Table S6). Because the average contribution rate of Bio11, Bio18, Bio7, and elevation to the Maxent model was higher than 10% and the accumulated contribution rate of these four variables was >80%, these four variables were classified as dominant variables affecting the distribution of sour jujube in different periods. The present potential distribution of sour jujube predicted in Maxent was highly consistent with the actual species' distribution (Wu, 1982). During the LGM, the range of sour jujube contracted considerably and showed a southward range shift. Peninsula) and around Bohai gulf. In the future (2050s and 2070s), the area of the total suitable area of sour jujube decreased but that of the higher suitable area increased. The distribution range of sour jujube clearly showed a northward shift to high-latitude areas in future environmental conditions ( Figure 5 and Table S6).

| Nucleotide variation of sour jujube in China
With the development of sequencing and analytical technology, genome data, such as genome-resequencing and RAD-seq (Restriction site-associated DNA sequencing) data containing tens of thousands of SNPs have been widely applied in population genetic research (Lange et al., 2021;Mu et al., 2020). Using genome-resequencing SNP data, the nucleotide variation of sour jujube was found to be higher than that of Chinese jujube (Guo et al., 2021;Huang et al., 2016). Though the number of SNPs in the present study (237 SNPs in five single-copy nuclear gene markers) was small, but the basic population genetic information of sour jujube was revealed to some extent and added to our acknowledgement of this important plant species. Future genome data may apply in the genetic-related study of sour jujube, such as the domestication history of sour jujube and the underlying mechanism for metabolite differences between sour jujube and Chinese jujube.
A series of factors may contribute to the nucleotide variation of plant species, including representative sampling, natural selection, mating system, and demographic history (Wright & Gaut, 2005). Representative sampling can be ruled out to a great extent because the sampling area of this study covered the main distribution range of sour jujube in China, especially the areas from the Loess Plateau to the Taihang Mountains. Natural selection can also be excluded because the test statistics indicated no significant deviation from neutrality. Sour jujube was obligated to outcross and inbreeding seldom occurred in natural populations (Zhang, 2013), which contributed to a relatively high recombination rate and effective population size, and ultimately to a high level of nucleotide variation (Charlesworth, 2003). Considering how living organisms responded to global climate fluctuations, specially the glacial-interglacial cycles in the Quaternary (Qiu et al., 2011), we primarily speculated that sour jujube may undergo a process of population contraction and expansion during the evolutionary history. But the BSP analysis revealed that the effective population size of sour jujube underwent a continuous increase.
This paradox can be clarified that the Quaternary glaciation had little effect on areas in Central China, especially the Lose Plateau (Zheng et al., 1998). In other central distribution areas of sour jujube, such as the Taihang Mountains, the population of plant species often migrated to high-altitude areas to prevent severe bottlenecks from occurring, and maintained a higher effective population size (Hewitt, 2000). Consequently, the mating system and complex population dynamics have acted to retain a high level of genetic diversity in sour jujube distributed in China.

| Genetic structure of sour jujube in China
In the present study, AMOVA analysis revealed that most of the total variation was attributed to variation within populations and a high level of genetic differentiation was detected (F st = 0.197, Table 3). This indicated that gene flow among populations was restricted to a relatively low level. Zhang et al. (2013) investigated the genetic structure of 3 sour jujube populations along the Yellow River and found a high level of gene flow and a low level of genetic differentiation among these populations. They speculated that seed dispersal depending on the water flow of the Yellow River may accelerate gene flow among these populations.
They further investigated the genetic structure of 34 sour jujube populations and also found moderate differentiation (F st = 0.091; Zhang et al., 2015). The populations sampled in Zhang et al. (2015) were mainly distributed along the reach of some rivers, such as the Yellow River, the Jing River and the Luo River, seed dispersal by the river over long distance was the most probable explanation for gene exchange. However, there was no obvious connection among these populations sampled in the present study. Furthermore, the Mantel test revealed a significant correlation between geographic, environmental distance, and genetic distance. Differences in the environmental conditions and long distance resulted in a relatively high level of genetic differentiation among the sampled populations in this study. The long-distance dispersal capability and the vitality of sour jujube pollen were relatively low (Shao et al., 2020). In general, it appears that forces such as isolation and adaptation that tend to increase genetic differentiation have been much stronger than homogenizing forces such as gene flow (Savolainen et al., 2007).  Zhao et al. (2021) compared 48 parameter combinations in R package ENMeval and Biomod2 to predict the distribution demography of sour jujube and found the Maxent model was optimal. In the present study, Maxent analysis yielded relatively accurate results that could be utilized to interpret the population demography of sour jujube corresponding to various environmental conditions in different periods. It was found that Bio11 contributed the most to the Maxent modeling, which reflected that low temperature was the main factor limiting the distribution of sour jujube. Sour jujube was not sensitive to soil conditions and its tolerance to barren soil was high Zhao, 2016), which resulted in the lowest contribution of soil to Maxent modeling. The effects of Quaternary glaciation-interglaciation cycles on the distribution pattern of plant organisms in the Northern Hemisphere have been substantially documented (Hewitt, 2004;Qiu et al., 2011). It has been proposed that temperate forests retreated southward to approximate 30°N during the Quaternary glacial periods based on paleovegetation data from East Asia and the populations located in northern areas must have recolonized from southern glacial refugia (Qian & Ricklefs, 2001 (Zheng et al., 1998). The mountainous areas were affected little by the Q uaternary glacier and the climatic condition was relatively stable in the high altitudinal area. Therefore, appropriate climate conditions in these regions can serve as the highly suitable areas for sour jujube to survive the glaciation. During MH and LIG, the distribution of sour jujube showed a northward recolonization and expansion resulting from the improvement of climatic conditions (Yu et al., 2000;Zhao & Piperno, 2000). In the future (2050s and 2070s under 3 RCP conditions), the distribution of sour jujube in China showed northward migration to high latitudinal areas.

| Population demography of sour jujube in China
Compared with the present distribution, the area of the higher suitable area increased but that of the total suitable area decreased.
This was consistent with the research that under warming climate conditions, the distribution area of plant species could reduce and shift to high latitudinal areas (Ji et al., 2020;Walther et al., 2005).
By simulating and predicting the distribution range of sour jujube in different periods, the historical causes of its formation and its future population dynamic responding to climate changes were de-

| Phylogeographic history of sour jujube in China
The effects of Quaternary climatic oscillations on the phylogeographic structure of species in the mid-to high-latitude regions of Europe and North America (Emerson et al., 2010;Ren et al., 2017),  Islam & Simmons, 2006 (Liu, 2013) and X. sorbifolium (Zhu, 2016), resulting from appropriate temperature and adequate precipitation (Wang et al., 2004). In Maxent analysis, it was revealed that the higher suitable distribution area of sour jujube in China located at the Loess Plateau in different periods.
During field sampling, it was found that sour jujube was distributed continuously in this area. Furthermore, nucleotide variation of central populations distributed in the Loess Plateau was much higher than that of the marginal populations (Table S8). Network analysis found that the chloroplast haplotype located in the Loess Plateau was more ancient. Combining these findings, we can speculate that a potential glacial refuge for sour jujube may locate in areas between the Loess Plateau and the Qinling Mountains, probably in areas adjacent to the Zhongtiao Mountains. Another higher suitable distribution area of sour jujube in LGM is located at areas adjacent to the Tianmu Mountains, which was also thought to serve as the main glacial refuge for temperate plant species in China, such as Ostryopsis davidiana Decne. (Tian et al., 2009), Juglans mandshurica Maxim. (Bai et al., 2018) and Ginkgo biloba L. . However, in the present we did not collect samples from this area, maybe because the environmental conditions were not suitable for sour jujube, such as high moisture.
These mountainous areas were important for the glacial survival of sour jujube and for preserving most of the extant nucleotide variation. This provides additional evidence for the in situ survival of plant species in Central China during the Quaternary glaciation and local expansion during interglacial or postglacial, which was similar to the scenarios revealed in other regions (Wang, Abbott, et al., 2010;Wang, Ikeda, et al., 2010). The most important refugia for sour jujube appear to have been located at the Loess Plateau were not collected in the present study, but based on the above analysis, we may speculate that the populations distributed in these areas may also migrate or colonize from the potential refugia or central populations revealed in the present study. Future research may include samples from these areas to explicitly test our speculation.

| CON CLUS ION
This study provided new insights into the genetic diversity, population structure, and demography of Z. jujuba var. spinosa. The genetic diversity of sour jujube was relatively high. Results of AMOVA showed that most of the total variation was attributed to variation within populations and high genetic differentiation among populations was detected. Geographic distance and environmental difference contributed to the genetic differentiation among populations. Investigation (equal).

ACK N OWLED G M ENTS
We thanked Pro. Wang Shengji and Wang Zhiling in College of Forestry, Shanxi Agricultural University in revising the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data presented in this study are openly available in NCBI GenBank with accession number MW371295-MW372934 for single-copy nuclear gene markers and NCBI SRA with accession number ON611607-ON611627 for chloroplast genome data.